Spontaneous development of 3-D structure in sheared granular flows 
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Computer simulations of sheared granular fluids, modeled as inelastic hard spheres, are presented 
which show signs of a uniquely three-dimensional instabilty. In the stable regime, a linear velocity 
profile, v x = ay, with shear rate a is established using Lees- Edwards boundary conditions. In 
the unstable regime, the velocity profile aquires a dependence on the third dimension of the form 
v x = ay + sin(2nz/L) in a cubic box with sides of length L. An analysis of the linearized Navier- 
Stokes equations shows the presence of an instability and gives a simple expression for the critical 
wavevector which is quantitatively consistent with the results of simulations and which indicates 
that the instability persists at low densities. 

PACS numbers: 45.70.Mg,47.20.Ft,83.10.Rs, 83.50.Ax 
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One of the intriguing aspects of fluidized granular flows 
is that they exhibit a rich phenomenology of clustering, 
segregation and pattern formation under different flow 
conditions (see e.g. |jj , [2| , |3j ) . Understanding the mecha- 
nisms that give rise to this behavior is one of the principle 
challenges in the field of granular physics. The minimal 
microscopic model of a granular fluid as hard spheres 
which lose energy upon collision is enough to give rise 
to many effects mentioned above. For example, com- 
puter simulations using this model of an isolated fluid 
shows that it is subject to both linear instabilities, which 
drive the formation of vortices 0,0,0, and nonlinear in- 
stabilities which lead to the formation of clusters 
Similar phenomena are observed in both simulation and 
experiment of the more practically relevant sheared gran- 
ular flows, in which the flow speed varies in some di- 
rection perpendicular to the direction of flow. Most at- 
tention has been focussed on 2-D flows due to the rel- 
ative computational efficiency with which they can be 
simulated and ease with which they can be studied in 
experiments and the results visualized although a few 
experiments have probed the 3-dimensional structure in 
annular geometries It has long been known that 2-D 
flows show the formation of both large scale structures 
(shear bands) and small-scale ordering of the grains[lfjj. 
The analysis of the stability of shear flow is complicated 
by the convective terms which give rise to a coupling 
of hydrodynamic modes 11] for wavevectors that sam- 
ple the direction of the flow but it has been shown that 
2-D unbounded shear flows are linearly stable though 
subject to nonlinear instabilities that drive the pattern 
formation[10j. Bounded 2-D shear flows, with moving 
walls perpendicular to the velocity gradient and periodic 
boundaries in the direction of the flow, have been shown 
to be linearly unstable 0,0. Dense granular flows in 3 
dimensions show solid-like ordering into regular arra ys |9| 
reminiscent of that seen in elastic sheared fluids|14 | .|l5 | . 

In this work, I present the results of computer simula- 



tions of the minimal granular model subject to shear flow 
in an unbounded geometry using Lees-Edwards bound- 
ary conditions and show that the system develops 3- 
dimensional structure even at low densities where there 
is no sign of the solid-like ordering. In particular, the 
system develops an intermittant variation of the pecu- 
liar velocity (relative to the linear profile) as a function 
of position in the direction perpendicular to both the 
flow and the velocity gradient with clustering near the 
extrema of this profile. A linear stability analysis of the 
relevant hydrodynamic equations shows the presence of a 
linear instability towards perturbations with wavevectors 
in the third dimension and a simple stability criterion is 
derived. Interestingly, this instability appears to be very 
similar to one which occurs in the stability analysis of 
the HCS|L6|L_|l 7j. 6] but which is hidden by the nonlinear 
instabilities [8|. Based on these results it appears that 
3-D sheared granular flows are linearly unstable for all 
densities. 

The simulations described here of mechanically iden- 
tical atoms were performed using an event-driven algo- 
rithm in which atoms stream freely between collisions. 
Two atoms with positions and velocities q i, v i and 



collide when their separation 



a where a is the hard-sphere diameter. The relative ve- 
locity 

collision by 



v '„■ after the collision is related to that before the 
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where %j = ~c?ij/qij is a unit vector and the total mo- 
mentum is unchanged. The parameter a is the coefficient 
of (normal) restitution: elastic hard spheres correspond 
to a = 1 while for smaller values of a, energy is lost 
and the collisions are inelastic. More complex models al- 
low for velocity-dependent coefficients of restitution and 
non-conservation of energy related to the tangential ve- 
locities as well but the focus here is on demonstrating 
the phenomena in the simplest model. The simulations 
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FIG. 1: Temperature as a function of the number of collisions 
for a — 0.7 with 13500 grains, upper black line, 500 grains, 
lower black line, and a = 0.9, N = 13500. white line. The 
arrows mark the times at which snapshots were taken. The 
arrows, from left to right, mark the snapshots (a)-(d) of the 
other figures. 



may be incorrect for the larger system. When a is set 
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FIG. 2: Snapshots of the flow velocity as a function of posi- 
tion along the y axis. The velocities are scaled to the initial 
thermal velocity. Each point is one of the 13500 grains. 



were performed in a cubic cell with walls of length L 
using Lees-Edwards boundary conditions 18| which are 
periodic boundaries in the Lagrangian frame defined by 
the velocity field ~v (~r*) = ayx, where a is the shear 
rate. It is easy to show that the balance laws admit of 
a stationary, spatially homogeneous state with the col- 
lisional cooling balancing the viscous heating to give a 
constant temperature. This balance relates the temper- 
ature, applied shear rate and coefficient of restitution so 
that there are effectively only three independent dimen- 
sionless parameters which will here be taken to be the 
number of atoms, N, the size of the cell, L/<r, and the 
coefficient of restitution a. The temperature can always 
be set to ksT = 1 by adjusting the time scale and the 
shear rate is then fixed by the energy balance relation. 
The procedure in all simulations was to first equilibrate 
an equilibrium (a = 1 and so a — 0) liquid of the de- 
sired number of atoms and cell size by starting with a 
face-centered cubic configuration and allowing it to melt 
over a run of 5 x 10 7 collisions. Long before the end of 
these simulations, the system exhibits the expected pres- 
sure and diffusion constant of a hard sphere liquid. The 
shear rate and coefficient of restitution are then set to 
the desired values and the simulations continued. 

A simple global measure of the behavior of the system 
is provided by the temperature defined as the kinetic en- 
ergy calculated relative to the assumed linear shearing 
profile. Figure 1 shows the temperature as a function of 
time for a — 0.7. The temperature of a small system, 
made up of 500 grains, is constant whereas that of the 
larger system, containing 13500 grains, increases rapidly 
and then appears to oscillate. The large deviation from 
between the two is an indication that the assumed profile 



to 0.9, the larger system also appears to be stable. Fig- 
ure 2 shows a snapshot of the velocity in the x— direction 
as a function of position in the y-direction for systems 
over one cycle of temperature oscillation while Fig. 2 
shows the peculiar velocity v x — ay as a function of z 
for the same snapshots. The v x vs y plot shows a gener- 
ally linear profile but with clear signs of local structure 
which is presumably a manifestation of the kind of 2-D 
structure described previously jlfjj. The plot vs. z shows 
a clear sinusoidal variation with a time-varying ampli- 
tude that grows to many times the thermal velocity at 
its maximum. From the density of points in the figure, 
it is clear that there is some clustering at the positions 
of the maxima and minima of the curves. The snap- 
shots taken at the minimum of the temperature cycle 
show that the 2-D structure persists but that the 3-D 
structure is virtually absent. Examination of numerous 
snapshots confirms that this oscillatory behavior takes 
place throughout the duration of the simulation and is 
the source of the temperature oscillations. As would be 
expected from Fig. 1 , which shows that the temperature 
profiles for a small system of 500 grains at a = 0.7 and 
and a 13500 grain system with a = 0.9 are flat, neither 
of these systems exhibits this structure. Simulations of 
the larger system for a = 0.6 show temperature varia- 
tions that are more than 5 times greater than those for 
a = 0.7. 

To understand why some of the simulations would 
spontaneously develop the 3-D structure, it is natural 
to examine the linear stability of the system. The local 
hydrodynamic fields, the density n(r , i), temperature 
T(V,£)and the velocity U (r satisfy the exact bal- 
ance equations flflj . |2(ij] 
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FIG. 3: Snapshots of the peculiar velocity in the x-direction 
as a function of position along the z axis. The velocities are 
scaled to the initial thermal velocity. Each point is one of the 
13500 grains. The white line is a fit to a sin-function with 
wavevector fixed at 2n/L. 
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where P is the pressure tensor, Q is the heat flux vector 
and ^ < is the cooling rate. Introducing the peculiar 
velocity u (r ,t) = U (r,t) — a • V and using the 



fluxes and cooling rate derived from the Enskog equation 
by expanding about HCS[l^| 
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where p is the pressure, 77 is the shear viscosity, 7 is 
the bulk viscosity, n is the thermal conductivity and fi 
is a transport coefficient which goes to zero at a = 1 
and describes the fact that density gradients drive heat 
flow in the granular system. All of these quantities are 
functions of the local hydrodynamic fields. Expanding 
these fields about steady, uniform values as n ( r , t) = 
71q + eni ( r ,t) + where e is a fictitious small param- 
eter used to order a perturbative expansion of Eqs.J2EJ), 
one finds at order e° that the only nontrivial constraint 



is the energy balance equation 

ma = —Co (4) 

where t]q is the shear viscosity evaluated for the zeroth- 
order fields. This gives the constraint relating the shear 
rate to the cooling rate, and hence to the temperature 
and the coefficient of restitution. Note that this uniform 
solution is consistent with the periodic boundary con- 
ditions in the co-moving frame used in the simulations. 
Further details of the calculation, which is straightfor- 
ward, will be given elsewhere and here, only the main 
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steps and results summarized. At linear order, it is con- 
venient to switch to a Fourier representation but the sta- 
bility analysis remains nontrivial due to the linear mode- 
coupling which occurs via the gradient term arising from 
the convective terms in eq.{(3J ■ However, to investigate 
stability in the z— direction, we can begin by restricting 
attention to wavevectors in that direction k = kz in 
which case, several simplifications occur. First, the gra- 
dient term does not contribute, so that temporal stability 
is simply determined by the signs of the eigenvalues of the 
propagation matrix, the elements of which are quadratic 
functions of the wavevector. Second, the x— and y— 
components of the velocity, representing shear modes, 
decouple form the other variables and are purely decaying 
with time constant equal to the kinematic shear viscosity 

I 



where the shear rate has been eliminated using the en- 
ergy balance equation. Here, a non-numerical subscript 
indicates a derivative evaluated in the homogenous state, 
e.g. v n — {j^y) _ ■ Recall that /io vanishes for a = 1 
so, for small deviations from elastic hard spheres, the co- 
efficient on the right will be positive and, since ^q '' > 0, 
the critical wavevector exists if and only if the term in 
brackets is positive. It is easy to show by differentiat- 
ing the secular equation that the slope of the real part 
of the eigenvalue at the critical wavevector is positive so 
that the mode is unstable for wavevectors smaller than 
the critical wavevector. The remaining two possibilities 
for critical wavevectors correspond to complex-conjugate 
propagating modes which become unstable but they are 
only relevant at extremely small values of a and will not 
be considered further here. Evaluation of this expression 
for the parameters of the simulation using the transport 
coefficients of ref.0] gives k c a — 0.238. The instability 
will not be present for small systems since the smallest 
nonzero wavevector that can be sampled is 2tt / L so that 
this defines a maximum stable system size or minimum 
number of grains in a given geometry. For a cubic system, 
this gives a stability limit of N c — 4612 which is indeed 
between the sizes of systems examined here. Similarly, 
for fixed N — 13500, the instability exists for a < 0.804 
which is in agreement with the observations above. These 
calculations also show that the left hand side is positive 
for all a at zero density. Nevertheless, this calculation 
cannot be taken too seriously as the dependence of the 
transport properties on the shear rate has not been taken 
into account, although it might be hoped that the vari- 
ous ratios appearing in eq.© would be less sensitive to 



is vq — ^^Vo ■ The remaining three eigenvalues are the 
roots of a cubic secular equation with fc-dependant coef- 
ficients. Setting k = shows that two modes go to zero 
and the third is stable with a finite decay constant. The 
latter assures that small, homogeneous violations of the 
energy balance relation eq. 10} decay. The stability of the 
modes can be determined by looking for solutions of the 
secular equation for which the real part of the eigenvalue, 
A, is zero, or in other words A = for real £. Then, the 
secular equation splits into two equations, the real and 
imaginary terms, which serve to determine both £ and 
the critical wavevector k c for which the mode does not 
decay. As expected, there are three possible solutions. In 
the first, £ = and the critical wavevector k c satisfies 



3™0 t (0) fPn \ , 2 / k n 

= ( — K o - Mo J K (5) 
I 

this omission than the absolute values of the quantities. 

The expression for the critical wavevector given in 
eq.© is closely related to a similar instability that oc- 
curs in the hydrodynamic equations linearized about the 
HCS Hi|, [l2|,|3. Obviously, in the HCS there is no vis- 
cous heating to counteract the cooling of the gas, so the 
system is not in a steady state and hydrodynamic modes 
in the usual sense do not exist. However, the system can 
be mapped to a steady state by introducing a new time 
variable related nonlinearly to the laboratory time. This 
introduces source terms into the hydrodynamic equations 
which replace the two terms in eq. (5j involving the shear 
viscosity but the analysis is otherwise identical. As dis- 
cussed in detail by Brey et al(3, this linear instability is 
not seen in the HCS because it is always slower to de- 
velop than the shear instability. The shear instability 
does not exist in sheared granular fluids and so it ap- 
pears that the secondary instability identified in HCS is 
actually manifested. 

These results only give information about the stabil- 
ity of the uniform velocity profile. Once the instability 
begins to develop, the nonlinear couplings in the hydro- 
dynamic equations will become important and will gov- 
ern the long time behavior. (Burnett terms and other 
higher order contributions to the fluxes neglected here 
will also become important.) In the present case, it is 
not clear whether these might exert a stabilizing influ- 
ence that counteracts the linear instability, or whether 
the oscillations are due to a frustration of the system 
by the cubic simulation cell and periodic boundaries. A 
careful analysis of the effect of the nonlinear terms fol- 
lowing along the lines of those performed for HCS[2l|.[8| 
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is clearly needed. 

Very recently, simulations of 3-D shear flows have 
shown the presence of density variations along the third 
dimension (perpendicular to both the flow and the veloc- 
ity gradient) |22j .The simulations were performed using 
hard walls at y — ±L/2 which move in the indirection 
thus creating a varying flow velocity between them and 
are characterized by the several collision parameters spec- 
ifying the normal and tangential coefficients of restitution 
(governing the amount of energy and angular momentum 
lost in the collisions), for both the atom- atom and atom- 
wall collisions and are also complicated by the inevitable 
presence of boundary layers near the walls. It seems likely 
that this is a manifestation of the same phenomena. 
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